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Abstract. We present quasi-stationary simulations of the two-dimensional 
contact process with quenched disorder included through the random dilution of a 
fraction of the lattice sites (these sites are not susceptible to infection) . Our results 
strongly indicate that the static exponents are independent of the immunization 
fraction. In addition, the critical moment ratios m = (p 2 )/(p) 2 deviate from the 
universal ratio m = 1.328, observed for the non-dilluted system, to smaller values 
due to rare favorable regions which dominate the statistics. 
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1. Introduction 

Absorbing-state phase transitions, i.e, transitions from a fluctuating phase to an 
absorbing (trapped) state, are related to several nonequilibrium critical phenomena 
[3 [21 [31 H] such as chemical catalysis [5], interface growth [5J, epidemic spreading 
[7J, etc. The study of such transitions in spatially extended systems has been 
experimenting an ongoing interest, strengthened by recent experimental confirmations 
of absorbing-state phase transitions in a liquid crystal system [8], and in a sheared 
colloidal suspension [9]. Notwithstanding a complete classification of their critical 
behavior is still missing, it has been conjectured [101 E] that models with a positive 
one-component order parameter, short-range interactions and deprived of additional 
symmetries or quenched disorder belong generally to the universality class of directed 
percolation (DP), which is considered the most robust universality class of the 
absorbing-state phase transitions. 

The contact process (CP) [12] is one of the simplest and most studied models of 
the DP universality class. Of particular interest is how spatially quenched disorder 
affects its critical behavior [T21[T3]. Quenched disorder, in the form of impurities and 
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defects, plays an important role in real systems, and may be responsible for the rarity 
of experimental realizations, in spite of the ubiquity of the DP class [T5] . 

The so-called Harris' criterion [T7| states that quenched disorder is a relevant 
perturbation, from the field-theoretical point of view, if 



where d is the dimensionality and v±_ is the correlation length exponent of the 
pure model (In DP this inequality is satisfied in all dimensions d < 4, since v±_ = 
1.096854(4), 0.734(4) and 0.581(5), for d=l,2 and 3, respectively Q21 EE ID] ) . The 
first numerical studies of the CP with quenched disorder, introduced by the means 
of a random deletion of sites [2TJ [22] or bonds [23] (dilution) or by random random 
spatial variation of the control parameter [24] . confirmed that the disordered system 
does not belong to the DP class [3T], and also revealed that the critical spreading is 
logarithmic, not a power law [22] . In the subcritical regime, a Griffits-phase, with 
critical dynamics dominated by nonuniversal power laws was also reported [22] 113] . 
However, the contact process in a Voronoi-Delaunay lattice, which has an intrinsic 
quenched disorder in the distribution connectivity, belongs to the DP class [25]. 

Hooyberghs et al. |26[ 127] . employed a strong-disorder renormalization group 
approach to conclude that the unusual critical behavior of the disordered system can 
be related to the random transverse-field Ising model, for sufficiently strong disorder. 
At such infinite-randomness fixed point, the scaling is activated, i.e, the temporal and 
spatial correlation lengths (£|| and respectively) are related by 



where ip is a universal exponent. For weak disorder they found nonuniversal critical 
exponents depending on the disorder strenght. 

More recently, Votja and Lee [28] used this activated dynamic scaling to show 
that the interplay between geometric criticality and dynamic fluctuations leads to a 
novel universality class, with the exponent ip equals to the fractal dimension of the 
critical percolation cluster of the diluted lattice. Previous numerical studies of the 
one-dimensional contact process with quenched spatial disorder, performed by Votja 
and Dickinson [29] also supported the activated exponential dynamical scaling at the 
critical point. Moreover, they found evidences that this critical behavior turns out 
to be universal, even for weak disorder. Novel strong disorder renormalization group 
calculations in a very recent paper by Hoyos 30] predict that the system is driven to 
the infinite-randomness fixed point, independently of the disorder strength. 

Thus, despite of a deeper understanding of the effects of quenched disorder 
in the critical contact process achieved in the last decade, a certain controversy 
remains: Do the static critical exponents change continuously with the degree of 
disorder J22J [13l [26] , or do they change abruptly to the values in the strong disorder 
limit corresponding to the universality class of the random transverse Ising model 
[29] [28] ? Obtaining the static exponents is a hard numerical task because at 
criticality the infinite disorder fixed point is characterized by a ultra-slow dynamics, 
p(t) ~ hx(t)~^' ( V J-^)^ leading to a unusual long relaxation towards the quasi-stationary 
(QS) values. Thus, in a tentative to shed some light into this issue, we present results 
of extensive large-scale QS simulations [31] of the two-dimensional diluted contact 
process. 

The balance of this paper is organized as follows: In the next section we review the 
definition of the contact process and describe the simulation method. In Sec. Ill we 
present our results and discussion, and Sec. IV is devoted to draw some conclusions. 



dv± < 2, 



(1) 




(2) 
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2. Model and methods 

The contact process is a stochastic interacting particle process in which the particles lie 
on the sites of a d-dimensional hypercubic lattice. Each site can be vacant or occupied. 
An empty site becomes occupied at a rate Xn/d, where n is the number of its occupied 
nearest-neighbors, while occupied sites become vacant at unitary rate [12j[T]. For a 
certain critical value A = X c {d), the model exhibits a continuous phase transition from 
an active state (with a positive density of sites occupied) to an absorbing configuration 
with all sites vacant, since none particle can be created from the vacuum. 

In the simulation we employ the usual scheme [T|: with probability p = A/(l + A), 
one nearest neighbor j of the selected site i is chosen at random and occupied if the 
site j is vacant. With complementary probability q = 1/(1 + A) the particle at site i 
is annihilated. At each step, the time is increased by At = l/N occ , where N occ is the 
total number of occupied sites. Moreover, occupied sites are sequentially selected at 
random from a list constantly updated, in order to improve efficiency. 

In the original diluted contact process (DCP) [3TJ a quenched disorder is 
introduced by labeling each site as diluted or nondiluted with probabilities x and 
1 — x, respectively. In the present work, we mark as diluted exactly a fraction T of the 
sites, in order to avoid undesirable extra fluctuations. Thus, DCP is the CP model 
restricted to the nondiluted sites since those diluted are never occupied. Notice that, 
for large systems, the dilution process used in this work is equivalent to those used in 
reference [22] . 

Stationary analysis nearby the critical point of systems with transitions to 
absorbing configurations are hard to be done due to very strong finite size effects. 
Indeed, the unique real stationary state of finite systems is the absorbing configuration. 
A common alternative to avoid this difficulty is to restrict the averages to the surviving 
samples and proceed a finite size analysis. This procedure [T], involves careful scrutiny 
in the data analysis which is not always free of ambiguities or misinterpretations [4]. 
In order to circumvent such difficulties, we employ a simulation method that yields 
quasistationary (QS) properties directly, the QS simulation method |31j . The method 
is based in maintaining, and gradually updating, a list of M configurations visited 
during the evolution; when a transition to the absorbing state is imminent the system 
is instead placed in one of the saved configurations. Otherwise the evolution is exactly 
that of a conventional simulation. By this procedure one obtains an unbiased sampling 
of the quasistationary distribution of the process. 

The simulations were performed as follows. Firstly, the list of configurations is 
incremented whenever the time increases by a unity up to a list with M = 1000 
configurations is achieved. Secondly, a configuration of the list is randomly chosen 
and replaced by the current one with a given probability p rep . We used a large value 
of p re p — 0.5 for a initial relaxation period, more precisely for t < L where L is 
the linear system size, in order to speed up the erasing of the memory of the initial 
condition. Also, p rep = 0.005 was adopted for the remaining of the simulation. 

3. Results and Discussion 

We studied lattice sizes varying from L = 20 to L = 640, averaging over 200 to 300 
independent realizations of disorder, of duration up to t = 2 x 10 8 . Averages were 
taken after a relaxation time t r — 10 8 . The larger sample sizes and longer run times 
apply to the larger L values. 
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The first step in analyzing our results is to determine, for each dilution value 
studied, the critical creation rate A(T). For this purpose we study the number of 
active particles, n(t) via spreading analysis. The critical value A c is then defined 
as the smallest A supporting asymptotic growth (vide Fig.l). This criterion avoids 
misinterpretations associated to the effects due to the Griffiths phase, in which power 
laws in n(t) are observed (22] . 

io 2 
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Figure 1. Spreading of activity from a single seed. Dilution rate: Y = 0.02. 

In Fig. 2 we show a log-log plot of the critical quasistationary density of active 
sites p, as function of L for dilutions ranging from 0.02 to 0.30. At criticality, 
such quantity decays as a power law, p ~ L~^l v ^. (This permits us to check the 
critical values A c obtained from the spreading analysis). The values for the critical 
exponent /3/v±, obtained from linear least-squares fits to the last four points of the 
data are shown in Table I. Our results suggest that the critical exponent ratio fi/v±_ 
is independent of the amount of dilution, with j3/v = 0.95(2), at least for dilutions 
r > 0.05. 

Table 1. Exponent ratio j3/v± for several dilution values. "^Present work using 
L < 640. 6 ) Taken from [22] where L € [8, 128] was used. 
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Figure 2. Quasistationary densities of active sites versus system size L. From 
top to bottom: nondiluted CP (dotted line), T = 0.02 (dashed), and T = 
0.05,0.10,0.20 and 0.30. 



Now we turn to the dynamical exponent z = v\\jv^. In the nondiluted CP, the 
lifetime of the QS state (which we take as the mean time between two attempts to 
absorbing in the QS simulation), follows r ~ L z . On the other hand, in the activated 
dynamical scenario, such power-law scaling is replaced by lnr ~ , with ip being 
an universal exponent. In other words, the critical exponent z is formally infinity in 
this scenario. Early works [221 126] found a nonuniversal power-law behavior, with an 
exponent z varying continuously in direction to the strong disorder values. Otherwise, 
our results, shown in Fig. 3, reveal that at criticality the lifetime behavior is not a 
power law, clearly diverging with an increasing slope for all T > 0.05. This suggest 
that the activated scaling emerge even for weak disorder, enforcing the universality 
hypothesis. Furthermore, applying the activated scaling to the the data for the last 
four points for highest dilution, T — 0.30 furnishes the value of tp ~ 0.48(7), consistent 
with the values of the exponent ip in the range 0.42 < t/j < 0.50 found in the literature 
for the random transverse Ising model 32, 33]. For weaker disorder, our results do not 
permit to distinguish between a scenario with continuously varying tp to a possible 
crossover to a universal value at larger system sizes, as predicted by the activated 
dynamics. 

Another consequence of the activated dynamics scenario is that the distribution 
of the observables become broader, implying that the averages are dominated by the 
rare events in which the process is locally supercritical. Thus, some quantities such 
as moment ratios of the order parameter (which converges to universal values in the 
nondiluted CP [34]) exhibit non-self-averaging properties, in the sense that they do not 
converge to limiting values even when L — > oo [2 2 j, I35| l36]. This is exemplified in Fig. 
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Figure 3. Quasistationary lifetime r versus system size L. Dilution rates: 
T = 0.05, 0.10, 0.20, 0.30, from bottom to top. Inset: dT/dL for T = 0.05. 



4 where the moment ratio m = (p 2 ) / (p) 2 is plotted as a function of the system size 
L. The effects of the rare 'favorable' regions dominate the statistics, and the moment 
ratio is drifted from the DP value of m ~ 1.328 to smaller values (in the favorable 
regions the process is locally supercritical, and m — > 1 in the limit A — > oo). We 
observe that for high dilutions the effects of the rare regions become observable even 
for modest system sizes, while for low dilution the effect only appear at considerable 
larger system sizes. Notice that for the dilution V = 0.02 these effects are not evident 
for the system sizes we used, what may justify the difference in the exponent ratio 
/3/v± for the smallest value of V. 

4. Conclusions 

We performed extensive large-scale simulations of the two-dimensional contact process 
with dilution. The dilution is known to change the critical behavior of the contact 
process. Our results indicate that the novel static exponents do not depend on 
the amount of dilution, and we present numerical evidences that the apparent 
nonuniversality observed in early works was due to finite-size effects. Our findings 
are in agreement with recent simulational results for the one-dimensional contact 
process with quenched disorder [28], and with strong disorder renormalization group 
results |30j . On the other hand, our results cannot exclude a nonuniversal variation 
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Figure 4. Moment ratio m for dilutions T = 0.30,0.20,0.10,0.05 and 0.02 from 
bottom to top. The dashed line represents the nondiluted V = 0.00 case. 



of the exponent ip with the disorder strength. Finally, the critical moment ratios 
m = (p 2 )/(/?} 2 deviate from the the universal ratio m = 1.328 of the non-diluted 
system to smaller values, due to rare favorable regions which dominate the statistics. 
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